/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm/common.hxx"
#include "uhm/util.hxx"
#include "uhm/mat.hxx"

#include "uhm/mat/linal.hxx"

namespace uhm {


  // Mat Ops
  // -------
  void random_matrix(Mat_Linal_ &A) {
    if (!A.is_flat_buffer_null()) {
      FLA_Random_matrix(~A.flat());
      //FLA_Set( FLA_ONE, ~A.flat() );
    }

  }
  void random_spd_matrix(Mat_Linal_ &A) {
    if (!A.is_flat_buffer_null()) {
      Int_ uplo = uplo2linal(A.get_uplo());
      FLA_Random_spd_matrix(uplo, ~A.flat());
      FLA_Hermitianize(uplo, ~A.flat());
    }
  }

  void norm1(Mat_Linal_ &A, Double_ &val) {
    if (!A.is_flat_buffer_null()) {
      FLA_Obj norm;
      FLA_Obj_create(type2linal_real(A.get_type()), 1, 1, 0, 0, &norm);
      FLA_Norm1(~A.flat(), norm);
      switch (A.get_type()) {
      case Mat_Base_::FLOAT  :
      case Mat_Base_::FCMPLX :
        val = *((float*)FLA_Obj_buffer_at_view( norm )); break;
      case Mat_Base_::DOUBLE :
      case Mat_Base_::DCMPLX :
        val = *((double*)FLA_Obj_buffer_at_view( norm )); break;
      }
      FLA_Obj_free(&norm);
    }
  }

  void norm_inf(Mat_Linal_ &A, Double_ &val) {
    if (!A.is_flat_buffer_null()) {
      FLA_Obj norm;
      FLA_Obj_create(type2linal_real(A.get_type()), 1, 1, 0, 0, &norm);
      FLA_Norm_inf(~A.flat(), norm);
      switch (A.get_type()) {
      case Mat_Base_::FLOAT  :
      case Mat_Base_::FCMPLX :
        val = *((float*)FLA_Obj_buffer_at_view( norm )); break;
      case Mat_Base_::DOUBLE :
      case Mat_Base_::DCMPLX :
        val = *((double*)FLA_Obj_buffer_at_view( norm )); break;
      }
      FLA_Obj_free(&norm);
    }
  }

  void copy(Mat_Linal_ &A, Mat_Linal_ &B) {
    if (!A.is_flat_buffer_null())
      FLA_Copy(~A.flat(), ~B.flat());
  }

  void copyt(Mat_Linal_ &A, Mat_Linal_ &B) {
    if (!A.is_flat_buffer_null())
      FLA_Copyt(FLA_TRANSPOSE, ~A.flat(), ~B.flat());
  }

  void axpy(Int_ alpha, Mat_Linal_ &A, Mat_Linal_ &B) {
    if (!A.is_flat_buffer_null())
      FLA_Axpy(int2fla(alpha), ~A.flat(), ~B.flat());
  }
  
  // UHM Ops
  // -------
  void random_matrix(UHM_Linal_ &U) {
#pragma unroll(UHM_MAT_LINAL_OPS_UNROLL)
    for (Int_ i=UHM_Base_::KIND_BT;i<(UHM_Base_::KIND_BB+1);++i) {
      UHM_Linal_::_B_ &B = U.B(i);
      if (!B.is_flat_buffer_null())
        random_matrix(B);
    }

#pragma unroll(UHM_MAT_LINAL_OPS_UNROLL)
    for (Int_ i=UHM_Base_::KIND_A_BEGIN;i<UHM_Base_::KIND_A_END;++i) {
      UHM_Linal_::_A_ &A = U.A(i);
      if (!A.is_flat_buffer_null())
        random_matrix(A);
    }

    //U.set_identity();
  }

  void random_spd_matrix(UHM_Linal_ &U) {
#pragma unroll(UHM_MAT_LINAL_OPS_UNROLL)
    for (Int_ i=UHM_Base_::KIND_BT;i<(UHM_Base_::KIND_BB+1);++i) {
      UHM_Linal_::_B_ &B = U.B(i);
      if (!B.is_flat_buffer_null())
        random_matrix(B);
    }

    // assume that this function is only used for cache
    if (!U.ATL().is_flat_buffer_null())
      random_spd_matrix(U.ATL());

    if (!U.ABR().is_flat_buffer_null())
      random_spd_matrix(U.ABR());
    
    //U.set_identity();
  }

  // BLAS 1 and 2 use flat matrices...
  void schur_x2r(UHM_Linal_ &U) {
    Int_ uplo = U.ABR().get_uplo();
    FLA_Gemm(FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE, FLA_ONE,
             ~U.ABR().flat(), ~U.XB().flat(), FLA_ZERO, ~U.RB().flat());
  }

  void copy(UHM_Linal_ &UA, UHM_Linal_ &UB) {
    copy_A(UA, UB);
    copy_B(UA, UB);
  }

  void copy_B(UHM_Linal_ &UA, UHM_Linal_ &UB) {
#pragma unroll
    for (Int_ i=UHM_Base_::KIND_B_BEGIN;i<UHM_Base_::KIND_B_END;++i) {
      UHM_Linal_::_B_ &S = UA.B(i);
      UHM_Linal_::_B_ &T = UB.B(i);
      if (!S.is_flat_buffer_null() && !T.is_flat_buffer_null())
        copy(S, T);
    }
  }

  void copy_A(UHM_Linal_ &UA, UHM_Linal_ &UB) {
#pragma unroll
    for (Int_ i=UHM_Base_::KIND_A_BEGIN;i<UHM_Base_::KIND_A_END;++i) {
      UHM_Linal_::_A_ &S = UA.A(i);
      UHM_Linal_::_A_ &T = UB.A(i);
      if (!S.is_flat_buffer_null() && !T.is_flat_buffer_null())
        copy(S, T);
    }

    copy(UA.P(), UB.P());
  }

  void axpy(Int_ alpha, UHM_Linal_ &UA, UHM_Linal_ &UB) {
    axpy_A(alpha, UA, UB);
    axpy_B(alpha, UA, UB);
  }

  void axpy_A(Int_ alpha, UHM_Linal_ &UA, UHM_Linal_ &UB) {
#pragma unroll
    for (Int_ i=UHM_Base_::KIND_A_BEGIN;i<UHM_Base_::KIND_A_END;++i) {
      UHM_Linal_::_A_ &S = UA.A(i);
      UHM_Linal_::_A_ &T = UB.A(i);
      if (!S.is_flat_buffer_null() && !T.is_flat_buffer_null())
        axpy(alpha, S, T);
    }
  }

  void axpy_B(Int_ alpha, UHM_Linal_ &UA, UHM_Linal_ &UB) {
#pragma unroll
    for (Int_ i=UHM_Base_::KIND_B_BEGIN;i<UHM_Base_::KIND_B_END;++i) {
      UHM_Linal_::_B_ &S = UA.B(i);
      UHM_Linal_::_B_ &T = UB.B(i);
      if (!S.is_flat_buffer_null() && !T.is_flat_buffer_null())
        axpy(alpha, S, T);
    }
  }


  // LINAL Interface
  // ---------------

  // diag is scalar operation
  Int_ apply_diag_scal ( Int_ side, Int_ conj, Int_ uplo,
                         Mat_Linal_ &A ) {
    return linal::host::internal::apply_diag_scal( side2linal(side), conj2linal(conj),
                                                   uplo2linal(uplo), 
                                                   ~A.flat() );
  }
  Int_ apply_diag_scal ( Int_ side, Int_ conj,
                         Mat_Linal_ &A, Mat_Linal_ &B ) { 
    return linal::host::internal::apply_diag_scal( side2linal(side), conj2linal(conj),
                                                   ~A.flat(), ~B.flat() );
  }
  Int_ apply_diag_inv_scal ( Int_ side, Int_ conj,
                             Mat_Linal_ &A, Mat_Linal_ &B ) {

    return linal::host::internal::apply_diag_inv_scal( side2linal(side), 
                                                       conj2linal(conj),
                                                       ~A.flat(), ~B.flat() );
  }

  Int_ apply_pivots( Int_ side, Int_ trans, Mat_Linal_ &p, Mat_Linal_ &A ) {
    return linal::apply_pivots( side2linal(side), trans2linal(trans),
                                ~p, ~A );
  }


  Int_ chol_nopiv  ( Int_ uplo, Mat_Linal_ &A ) {
    return linal::chol_nopiv( uplo2linal(uplo), ~A );
  }
  Int_ chol_piv    ( Int_ uplo, Mat_Linal_ &A, Mat_Linal_ &p ) {
    return linal::host::internal::chol_piv( uplo2linal(uplo), ~A.flat(), ~p.flat() );
  }
  Int_ chol_incpiv ( Int_ uplo, Mat_Linal_ &A, Mat_Linal_ &p ) {
    return linal::chol_incpiv( uplo2linal(uplo), ~A, ~p );
  }


  Int_ ldl_nopiv   ( Int_ uplo, Mat_Linal_ &A ) {
    return linal::ldl_nopiv( uplo2linal(uplo), ~A );
  }
  Int_ ldl_piv     ( Int_ uplo, Mat_Linal_ &A, Mat_Linal_ &p ) {
    return linal::host::internal::ldl_piv( uplo2linal(uplo), ~A.flat(), ~p.flat() );
  }
  Int_ ldl_incpiv  ( Int_ uplo, Mat_Linal_ &A, Mat_Linal_ &p ) {
    return linal::ldl_incpiv( uplo2linal(uplo), ~A, ~p );
  }


  Int_ lu_nopiv    ( Mat_Linal_ &A ) {
    return linal::lu_nopiv( ~A );
  }
  Int_ lu_incpiv   ( Mat_Linal_ &A, Mat_Linal_ &p ) {
    return linal::lu_incpiv( ~A, ~p );
  }
  Int_ lu_piv      ( Mat_Linal_ &A, Mat_Linal_ &p ) {
    return linal::lu_piv( ~A, ~p );
  }



  Int_ trmm  ( Int_ side, Int_ uplo, Int_ trans, Int_ diag,
               Int_ alpha, Mat_Linal_ &A, Mat_Linal_ &B ) {
    return linal::trmm( side2linal(side), uplo2linal(uplo),
                        trans2linal(trans), diag2linal(diag),
                        int2fla(alpha), ~A, ~B );
  }
  Int_ trsm  ( Int_ side, Int_ uplo, Int_ trans,
               Int_ diag, Int_ alpha, Mat_Linal_ &A, Mat_Linal_ &B ) {
    return linal::trsm( side2linal(side), uplo2linal(uplo), trans2linal(trans),
                        diag2linal(diag), int2fla(alpha), ~A, ~B );

  }
  Int_ gemm  ( Int_ transa, Int_ transb, Int_ alpha,
               Mat_Linal_ &A, Mat_Linal_ &B, Int_ beta, Mat_Linal_ &C ) {
    return linal::gemm( trans2linal(transa), trans2linal(transb),
                        int2fla(alpha), ~A, ~B, int2fla(beta), ~C );
  }
  Int_ gemm_l( Int_ transa, Int_ transb, Int_ alpha,
               Mat_Linal_ &A, Mat_Linal_ &B, Int_ beta, Mat_Linal_ &C ) {
    return linal::gemm_l( trans2linal(transa), trans2linal(transb),
                          int2fla(alpha), ~A, ~B, int2fla(beta), ~C );
  }
  Int_ gemm_u( Int_ transa, Int_ transb, Int_ alpha,
               Mat_Linal_ &A, Mat_Linal_ &B, Int_ beta, Mat_Linal_ &C ) {
    return linal::gemm_u( trans2linal(transa), trans2linal(transb),
                          int2fla(alpha), ~A, ~B, int2fla(beta), ~C );
  }
  Int_ syrk  ( Int_ uplo, Int_ trans, Int_ alpha,
               Mat_Linal_ &A, Int_ beta,  Mat_Linal_ &C ) {
    return linal::syrk( uplo2linal(uplo), trans2linal(trans),
                        int2fla(alpha), ~A, int2fla(beta), ~C );
  }
  Int_ herk  ( Int_ uplo, Int_ trans, Int_ alpha,
               Mat_Linal_ &A, Int_ beta,  Mat_Linal_ &C ) {
    return linal::herk( uplo2linal(uplo), trans2linal(trans),
                        int2fla(alpha), ~A, int2fla(beta), ~C );
  }
}
